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ABSTRACT 


I 1 


A variety  of  simple  analytical  models  for  increasing, 
decreasing  and  "bath  tub"-type  failure  rates  are  discussed. 
The  purpose  of  this  thesis  is  to  develop  analytical  hazard 
representations  for  use  in  reliability  and  maintainability 
studies,  and  to  evaluate  them  in  use  for  data  analysis. 
Verification  of  the  model  was  accomplished  by  computer  simu- 
lation. They  were  applied  to  human  mortality  and  other 
failure  time  data. 
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I.  INTRODUCTION  AND  SUMMARY 


A.  INTRODUCTION 

The  failure  rate  function,  or  hazard  function  (hazard 
for  short)  may  be  described  as  the  conditional  probability 
of  an  equipment's  failing  at  operating  age  t,  having  sur- 
vived to  that  age.  The  reliabilities  of  a variety  of  elec- 
tronic and  mechanical  items  are  conveniently  and  naturally 
described  in  terms  of  the  appropriate  hazard  function,  and 
so  is  the  longevity  of  human  beings.  The  term  force  of 
mortality  replaces  hazard  in  the  latter  context. 

This  paper  is  devoted  to  a study  of  several  simple 
analytical  representations  for  hazard  functions.  These 
representations  are  in  turn  based  upon  representations  of 
random  variables  having  certain  required  properties,  in 
terms  of  others  having  familiar  distributions — in  particular 
the  exponential.  Similar  ideas  are  due  to  Tukey  [1],  and 
recently  have  been  examined  by  Parzen  [2].  The  hazard 
representations  proposed  are  quite  expeditiously  used  in 
simulation  studies,  e.g.  of  system  reliability  or  avail- 
ability in  terms  of  component  lifetimes.  They  may  also  be 
used  in  data  analysis  studies,  in  order  to  parsimoniously 
describe  data  sets  in  terms  of  perturbations  of  convenient 
and  familiar  standard  distributions.  Their  use  in  data 
analysis  and  simulation  is  also  described  in  Gaver,  Laven- 
berg,  and  Price  [3],  and  in  Gaver  and  Chu  [4], 
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B.  SYSTEM  FAILURE  PATTERNS 

It  is  plausible  to  think  that  the  time  series  of  fail- 
ures in  a system  may  involve  these  stages. 

Early  failures.  There  may  be  a relatively  large  number 
of  failures  soon  after  a system  is  introduced  because  of 
design  defects,  production  errors,  or  errors  stemming  from 
maintenance  personnel  inexperience.  This  situation  is 
characterized  by  a hazard  function  that  is  initially  large, 
but  that  decreases  with  time.  "Infant  mortality"  is  in 
evidence. 

Random  Failures.  Following  the  early  failure  period 
there  may  be  a period  during  which  failures  occur  at  an 
essentially  constant  rate  for  a rather  prolonged  time. 
During  this  period  the  hazard  function  is  nearly  constant, 
so  the  times  between  failures  are  close  to  being  exponen- 
tially distributed.  The  effect  of  age  or  wearout  is  not 
yet  apparent. 

Wearout  Failures.  Eventually  following  the  period 
during  which  a constant  hazard  is  evident  there  is  likely 
to  be  a period  of  ever-increasing  failure  rate  caused  by 
wearout  of  system  components. 

A graphical  representation  of  a hazard  function  that 
exhibits  the  behavior  described  is  given  below.  Note  that 
it  has  the  legendary  "bath  tub"  shape. 
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Fig.  1.  Bath  Tub  Shape  Curve 
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Some  comments  on  the  above  follow: 

The  term  "failure"  may  refer  to  an  event  that  is 
analogous  to  human  death,  after  which  the  entire  system 
is  replaced.  On  the  other  hand  repair  or  component  replace- 
ment  may  occur  after  failure:  the  system  is  only  repaired, 
not  entirely  replaced.  In  the  former  case,  a hazard  func- 
tion  of  the  kind  depicted  in  Figure  1 applies  to  each  system 
event  ("death");  when  the  system  is  installed  (or  is  born), 
that  hazard  operates  starting  from  scratch  at  t = 0 until 
system  failure  (human  death,  for  instance) , after  which  a 
similar  hazard  goes  into  effect,  starting  once  again  from 
zero.  In  the  latter  case,  in  which  repair  of  a component 
occurs,  a hazard  function  like  that  of  Figure  1 applies  at 
t = 0,  but  after  the  first  event  ("failure")  at  t^  a re- 
pair action  is  accomplished.  The  same  hazard  operates  for 
t ^ t^  until  the  next  event  at  t2  > t^,  and  so  on.  Inter- 
mediate situations  may  be  envisioned,  in  which  after  event 
n at  t the  hazard  governing  system  failure  n+1  starts 

at  t - t , 0 < t < t . 
n n n n 

Although  there  is  reason  to  assume  that  hazards  some- 
what like  that  of  Figure  1 occur  in  general  for  systems, 
the  possibility  exists  that  the  system  hazard  is  "bumpy" 
because  wearout  failures  of  components  or  subsystems  may 
well  occur  at  intermediate  times. 

If  the  theory  is  applied  to  systems  with  little  or  no 
wearout  propensity,  as  should  be  the  case  when  dealing  with 
computer  software  modules,  then  the  hazard  function  may 
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well  exhibit  the  initial  falloff  of  Figure  1 but  not  the 
rise  at  later  times.  In  fact,  a constant  decline  as  bugs 
are  found  and  removed  could  be  (optimistically)  anticipated 
for  software.  The  right-hand  side  of  the  bath  tub  vanishes, 
and  the  picture  is  that  of  a ski  slope. 


II.  ANALYTICAL  HAZARD  REPRESENTATIONS 

A.  MODELS  FOR  THE  HAZARD  FUNCTION 

In  this  section  mathematical  models  are  presented  for 
the  failure  rate  or  hazard  function.  Recall  that  the  hazard 
may  be  defined  as  follows. 

Definition.  Suppose  that  the  time  to  failure,  X,  is  a 
random  variable  with  distribution  function  F(x) , where 
F(0)  ■ 0;  the  latter  possesses  the  density  function  f(x), 
f ( x)  - dF/dx,  such  that  for  any  positive  x, 

x 

F (x)  = / f(y)  dy  . (2.1) 

0 


Then  the  hazard  function , 
given  by 

h (x) 


or  failure  rate  at  age 

f (x) 

’ 1 - F(x)  * 


x,  is 


(2.2) 


The  interpretation  of  h(x)dx  is  that  it  is  the  con- 
ditional probability  of  failure  in  the  interval  (x,  x + dx) , 
given  that  there  has  been  no  failure  up  to  age  x. 

Express  the  hazard  as 


h (x) 


dF/dx 

1 - FTxT 


h ( x)  dx  - ^ - f(xV  * ~ d{1°9[1  ~ F(x)]}  , 


it  then  follows  after  integration  that 
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(2.3) 


x 

F ( x)  - 1 - exp [-  / h(y)dy)  . 

0 


Thus  if  the  hazard  is  specified,  so  is  the  distribution 
function,  and  conversely. 

Note  that  if 

h(x)  - X > 0 , 0 < x 


then 


F(x)  - 1 - e“Ax 


0 < x , 


(2.4) 


so  a constant  hazard  function  implies  the  exponential  dis- 
tribution of  the  random  variable  X,  and  conversely. 

Obviously  a constant  hazard  representation  does  not 
describe  the  bath  tub  hazard  shape  of  Figure  1,  nor  does  it 
represent  a situation  in  which  hazards  decline,  possibly 
because  design  defects  or  "bugs"  are  occasionally  removed. 
Here  are  two  hazard  representations  likely  to  be  useful  for 
such  purposes . 


1.  A Bath  Tub  Model 

Define  the  random  variable  Z in  terms  of  X,  X 
being  exponent ially  distributed  with  mean  X-1,  as  follows: 


Z - 

G(X)  - 

XL  (X)  R (X) 

or 

m 

X*(X) » 

<t>  (X)  - L(  X)  R (X) 

where 

a) 

L(x) 

is 

concave 

X 

c 

L(0)  < 1,  L<~)  - 1, 

b) 

R(x) 

is 

convex 

in  x, 

R (0 ) - 1,  R (0 ) > R (°°)  . 
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Then  the  hazard  of  Z may  be  made  to  exhibit  a bath  tub 
shape,  as  in  Figure  1,  by  proper  choice  of  the  functions  L 
and  R. 


Example . Suppose 


L<*>  - a > o,  o < x 


(2.6) 


R(x)  * 


T"  Bx 


3 > 0 , 


Clearly, 


z ■ xL(x)R(x)  ■ 


ax 

+ ax 


ax 

+ ax 


x 

T7x 


is  a monotonically  increasing  function  of  x.  Furthermore, 
choose  a large  (e.g.  a « 10)  and  6 small  (e.g.  3 • 10  ^) 
Then  it  is  intuitively  clear  that  (i)  small  x-values 
transform  into  even  smaller  z-values,  e.g.  x - 1 corre- 
sponds to  z - 0.91  and  x » 2 corresponds  to  z ■ 1.90, 
but  (ii)  this  effect  dwindles  as  x increases,  so  x - 10 
corresponds  to  z * 9.8  and  x » 50  to  z-47.5  and  the 
z-values  closely  resemble  the  x's  percentage-wise,  but 
(iii)  as  x increases  still  further  the  z's  do  not  follow 
suit:  x ■ 10 3 corresponds  to  z ■ 500.  This  suggests 
that  if  x is  a value  assumed  by  X,  that  Z shares  the 
properties  of  X in  mid-range,  i.e.  for  intermediate 
x-values,  but  differs  from  X by  having  a disproportionate 
probability  of  assuming  small  values  (near  zero) , or  large 
values  (near,  but  less  than,  1/3).  Thus  the  hazard  of 
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Z will  appear  to  be  a "bath  tubbed"  version  of  X, 
particularly  if  X is  exponential. 


i 


We  focus  attention  on  the  representation  (2.6)  in  what 
follows,  mainly  for  analytical  and  computational  convenience. 
Of  course  there  are  many  other  possibilities,  such  as 

L(x)  - 1 - e”°x 

-fix  (2‘7) 

R(x)  - e **  , 

these  latter  may  be  adjusted  to  provide  sharper-edged  tubs 
than  can  (2.6),  but  iteration  of  (2.6)  may  be  induced  to 
accomplish  the  same  purpose. 

2 . A Decreasing  Failure  Rate  Model 

Define  the  random  variable  w in  terms  of  X,  X 
again  being  exponential  with  parameter  1 ^ : 

W - XT  (X)  (2.3) 

where  T(x)  is  an  increasing  function  of  x,  L(0)  ■ 1. 

Then  the  hazard  may  be  made  to  exhibit  a decreasing 
behavior. 


Example.  Suppose 


T (x) 

Then 

is  monotonic, 


- 1 + cx  , c > 0 , 0 <_  x . (2.9) 

z * x(l  + cx)  (2.10) 

and  small  x-values  lead  to  comparable  z- values 


i 

i 
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(especially  when  c is  small) , but  larger  x-values  are 
"amplified"  by  1 + cx  to  yield  increasingly  large  z-values. 

Attention  will  be  focused  upon  (2.9),  although  other 
possibilities  exist  that  accomplish  the  same  purpose,  namely 
that  of  lengthening  the  right  tail  of  the  distribution  of 
X (simulating  outliers,  for  instance)  while  leaving  the 
body  of  the  distribution  virtually  unchanged. 

B.  MATHEMATICAL  PROPERTIES  OF  THE  "BATH  TUB"  HAZARD  MODEL 

Various  analytical  properties  of  the  previously  described 
models  will  now  be  recorded.  These  provide  useful  insights 
into  the  behavior  of  the  random  variables  Z and  the  under- 
lying (generating)  variables  X. 

1.  Monotonicity;  Quantiles 

It  is  convenient  to  focus  on  monotonic  increasing 
trams  formations , i.e.  if 


z ■ G(z)  » zp  (z)  (2. 10) 

then  in  order  that  the  above  function  be  monotonically  in- 
creasing, dz/dx  > 0.  Observe  that  logarithmic  differen- 
tiation of  (2.10)  provides 


dz 

z 


<*'  (x) 

pjzr 


dx 


and  thus  dz/dx  >0  if  and  only  if 


1 

x 


+ 


: (x 

17 


> o 


(2.11) 


(2.12) 
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Alternatively,  the  condition  is,  in  terms  of  L(x)  and 
R(x) , 


1 . L*  (x)  . R'  (x)  v « 

x + “TTxT  + > 0 


(2.13) 


It  is  easily  seen  that  the  important  example  (2.6), 


*(x)  “ r+x_ax  * r~rjz  • 


yields  a monotonic  relationship  between  z and  x.  The 
fact  that  this  transformation  can  be  easily  and  explicitly 
inverted  (solved  for  x in  terms  of  z)  will  be  exploited 
subsequently. 

Of  course  if  z(x)  is  monotonically  increasing  then  so 
is  x(z)  , the  inverse  f met  ion.  The  events  (Z  £ z)  and 
(X  <_  x(z))  are  equivalent,  and  so 


P{Z  £ z}  - P{X  <_  x ( z ) } , (2.14) 

from  which  it  follows  that  if  x = x(p)  is  the  p*100% 

P 

quantile  of  X,  i.e. 

P{X  £ x ( p) } » p , (2.15) 

then 

P{Z  <_  z (p)  } * P{Z  <_  z(x(p))}  » p (2.16) 

and  so  z(p),  the  p*100%  quantile  of  Z is  simply  obtained 
from 

z(p)  ■ x (p)  4>  (x (p) ) = x ( p)  L(x (p) ) R(x(p) ) (2.17) 

In  other  words  we  very  easily  translate  from  (points  on) 
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the  Inverse  distribution  of  X to  the  inverse  distribution 
of  Z.  Explicit  representation  of  the  distribution  of  Z 
is  however,  not  often  easily  possible. 

2 . Hazard  and  Density  Function  Relationships 

In  order  to  investigate  the  relationship  between 
the  hazards  of  Z and  X,  begin  by  writing 

x(p) 

p ■ FY (x (p) ) - 1 - exp [-  / h (u) du]  (2.18) 

A q X 

or 

*(P> 

/ h (u)du  - - in(l-p) 

0 x 

Now  differentiate  with  respect  to  p to  find 

hx(x(p))  d^P)~  * I=p  (2.19) 

or 

hx(x(p))  ■ aifpT  ' i=p  - £x(x<p"  • it?  • «•«> 

here  hx  and  fx  are  the  hazard  and  density  functions  of 
the  r.v.  X.  The  relationship  (2.20)  holds  for  any  distri- 
bution, of  course. 

Differentiation  of  (2.5)  reveals  the  connection  between 
h_  and  h . From  (2.11) 

Z X 

- (<Mx(p))  + x(p)  ♦•(x(p))] 


dx(p) 

“^P 


(2.21) 
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Prom  (4.11),  applied  now  to  the  z -hazard,  there  results 

h, Utpl)  * h U(pl)  [»<»(P)>*»(P>  ♦•<*(?»)  . (2-22) 


SO 


h*(*(P>)  * hx(x(p))  <j)(x(pl  )"+"xfp)  -0'Tx'(p)7  (2*23) 


Multiplication  of  both  sides  by  1-p  then  shows,  in  view 
of  (4.11) , that  the  density  functions  are  similarly  related: 


fz(z(p) ) * fx(x(p) } *(x(p))"+“x7p)  Vlx(p)T  * 


Example . 

X is  exponentialU ) . Then 


hz(x(p))  = <Mx(P>)  + xTpT  ?T7xTp7T  (2*24) 


Now  use  the  specific  0(x)  of  (2.6)  : 


A,..»  ax  1 

*(x)  “ l'~+~  ax  * 1"+  Fx  ' 


or,  in  terms  of  logarithms, 


In  ( x)  = Jin  ax  - Jln(l  + ax)  - Jln(l  + Sx)  , 


so 


O'  (x)  _ 1 

0^7  * X 


a 6 1 ~ «8x ,, 

1 + ax  1 + £x  “ x(l  + ax) (1  + 6x) ' 


and 


0 ( X)  + X0 ' (x)  * 0 (x) 


1 2 + (a  + 6 

I (1  + ax)  (I  T 


B)x 

TxT 


(2.26) 


finally 
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hzU(p)) 


Ml  + ax(p))2  (1  + Bx(p))2 
oix(p)  T2  + ("of + "8 r-x (pT r 


(2.27) 


Although  this  expression  is  not  quite  explicit,  qualitative 
properties  of  hz  can  be  deduced  from  it. 

If  p * 0,  x(p)  ■ - j in(l-p)  + 0,  and  hence 


hz(z(p)) 


icx  x(p)  ' 


(2.28) 


lim  x (p)  h (z (p) ) ■ k— 

p > 0 


(2.29) 


Since  for  p 0 , 


and  hence 


z (p)  - ax  (p) 


x(p)  ~ [z (p) /a ] 


(2.30) 


there  results 


h (z  (p) ) 


2(z(p)  /a 


lim  /iTpT  h (z(p))  « — 
p -►  0 z 2 r a 


(2.31) 


This  shows  that  h (z)  t ~ as  z ♦ 0 , creating  the  left- 

z 

hand  end  of  the  bath  tub  of  Figure  1. 

If  p > 1,  x(p)  t »,  and  z(p)  t 1/6  so 


~ * ll  + af1 


1 


(2.32) 


lim  i—y  h (Z(p))  -A-a£ 

p-«  [ x ( p)  ]2  2 a + 


For  p * 1 


so 


and  thus 


1 " 6z<p)  ~ SiKxfpy 


a + 8 1 

x(p)  — zr  r^TiTFr 

h (2(P))  L__ 

2 ' a ' (1  - 8z) 


(2.33) 


(2.34) 


(2.35) 


Once  again  it  appears  that  the  hazard  rises  rapidly,  this 
time  as  x(p)  t « and  z(p)  t 8 S the  other  end  of  the 
bath  tub  is  thus  fashioned. 

If  p * 1 - e then  x(p)  * Then 


h,U(l  - e'1))  - ,(1  + VM2  [1  * a/XI2  (2-36) 

o/A  [2  + (a  + 8)/X] 


Tho  bath  tub  effect  is  presumably  achieved  by  choosing  a 
large  and  6 small.  Let  a -*■  0 and  8 0 independently 

in  (2.36);  it  is  clear  that  the  limiting  value  of  the 
hazard  is  X.  This  indicates  that  the  hazard  is  (approxi- 
mately) X for  middling  values  of  z. 


3.  An  Explicit  Formula  for  a Hazard 

The  expression  (2.6)  leads  to  the  relationship 

2 

2(p)  * TT  — Sx{pf ) trVsxTp) T ' 


(2.37) 
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and  the  latter  may  be  explicitly  inverted  by  solving  a 
quadratic  equation.  The  result  is 


x (p) 


(g+B)z(p)  + Aa+ B)2  zS(p)  + 4z(p)  ot  (1  - gz  (p) ) 
2q(l_-8z  (p) } 


/ O O O \ 


Now  a direct  differentiation  of  this  expression  and  invo- 
cation of  (2.20)  produces  the  expression 


h,<*> 


2a'(l-ff*T 


8{(g+8)z  + >/(q+8)2zZ  + 4qz(l-gz)  + (a  + g) 
1 — gz 


{ (g-g)  2z  + 2<x}J(ol+&)2z2  + 4gz  ( 1-gz) 
(g  + g}^  + 4a* (1  - gz) 


(2.39) 


This  form,  while  explicit,  provides  no  particularly  useful 

insights;  the  bath  tub  end  shapes  already  noted  in  (2.31) 

and  (2.35)  can  be  deduced  directly  from  (2.39). 

Some  graphical  plots  of  h are  presented  below.  They 

z 

illustrate  the  behavior  of  the  present  hazard  representation 
in  a more  understandable  fashion  than  does  the  formula 
itself. 

4.  An  Explicit  Formula  for  the  Failure  Time  Distribution 
Because  x and  z are  monotonically  related 
through  (2.37)  we  have 
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P{Z  £ z } 

* p | x £ 


♦ '/(a+B)  2z2  4az  ( l-$z)  I 

I 


„ _ (a+B) z 

x imn^zT 


1 - exp  \ -X 


(a+8)z  + >/(a+B)2z2  + 4az(l-Bz) 
2a ( 1 - Bz) 


(2.40) 


Again  the  explicit  formula  seems  unproductive  of  insights. 


C.  MATHEMATICAL  PROPERTIES  OF  THE  DECREASING  FAILURE 
RATE  MODEL 

1.  The  Hazard  Behavior 

The  expression  (2.23)  c an  be  applied  to  deduce  the 
hazard  function  of  the  representation  (2.8),  advertised  to 
produce  a decreasing  failure  rate.  There  we  specified 

$(x)  = T(x)  * 1 + cx  , (2.41) 

and  thus,  from  (2.2  3) 

hw(w(p))  3 71  V CxTpiY"xfpTc  3 T + 2cxTpT  (2,42) 

Qualitative  properties  follow  easily. 

If  p -*■  0 , x(p)  + 0,  and 

hw(w(p))~X  . (2.43) 

Thus  the  hazard  is  approximately  X for  small  z. 

If  p * 1,  x(p)  f *,  and 

2 

w ( p)  c (x(p)  ] 

so 
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h (w (p)  ) ~ — • - — ' 

2 /cwTpT 


(2.45) 


which  clearly  decreases,  as  claimed.  It  may  be  inferred 
that  the  distribution  of  W appears  nearly  exponential, 
but  possesses  an  extraordinarily  long  right  tail — possibly 
the  result  of  outliers. 


Explicit  Formulas  for  the  Hazard  and  the 
Distribution  Function 

Direct  solution  of  the  quadratic  equation 


w * x(l  + cx)  * x + cx 


presents 


(2.46) 


which,  when  differentiated,  leads  to 


Vwl 


h (x) 
x 


(2.47) 


* 1 + 4cw 


The  distribution  function  is 


P{W  < w)  = P I X < x = 


r- 


1 - exp  -A 


I -> 


+ 4cw 


(2.48) 


This  distribution  bears  a close  family  resemblance  to  tAe 
Weibull  distribution  1 - F(w)  = exp{-k/w),  especially  for 
large  (right  tail)  values  of  w. 
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D.  AN  ALTERNATIVE  "BATHTUB"  HAZARD  REPRESENTATION 


The  simple  parametric  model  (2.6)  leading  to  a bathtub- 
shaped hazard  is  by  no  means  the  only  possibility.  We  next 
describe  another  simple  approach.  It  is  that  of  defining 
a hazard  function  having  an  appropriate  shape,  and  then 
deducing  the  corresponding  distribution  function,  and  a 
procedure  for  sampling  from  it,  rather  than  proceeding  in 
reverse  order,  as  before. 

Let  the  hazard  be  of  the  form 


h(z)  = g(z)  +•••+  k(z)  , (2.49) 

where  g(z)  >0  is  a decreasing  function  of  z such  that 
lim  g(z)  =0,  and  k(z)  is  an  increasing  function  of  z, 
such  that  (preferably)  k(0)  =0  and  k(°°)  = °°.  Such  a 
function  can  yield  a bathtubbed  hazard. 


Example . 

h (z)  = 2— + Bz  + A 


(2.50) 


A,  B,  a,  A all  positive. 

Clearly,  (2.50)  has  a generally  "bath  tub-like"  appear 
ance,  since 


if 


while  for 


h (z)  > 0. 


h'(z)  — B (2.51) 

(z  + a) ^ 

- | + B < 0,  then  h ' (0)  < 0 (2.52) 

z > zQ  = | - a (2.53) 


! 


... 
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Detailed  behavior  is  adjustable  by  choice  of  the 
parameters.  Now  the  distribution  function  of  time  to 
failure,  Z,  is  obtained  from  (2.50)  : 


i 


p{Z  > z}  « exp  j-  / h(x) dx  | 


-exp  - / (5TTJ 


A + Bx  + X)dx 


■ exp 


j - £ A Zn(l  + |)  + | z2  + Xz]  { 


n A . B 2,  -Xz 

<3-r7>  eJtp<'  i z e 


Fx(z)  F2(z)  F3(z)  , 


(2.54) 


where 


F^z) 


F2(z) 


F3(z) 


( )A 

\ a + z/ 

(-  I *2) 


exp 

-Xz 


■ e 


(2.55) 


All  of  the  above  are  recognized  as  being  the  complements  of 
distribution  functions.  In  effect,  the  distribution  of  Z 
is  that  of  the  minimum  of  three  independent  random  variables 


p{z  > Z } - P{X.  > z 


}-P{X,  > Z}*P{X.  > *}, 


(2.56) 
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having  the  distribution  ■ 1 - (i  - 1,2,3). 


This  fact  leads  directly  to  an  easy  procedure  for  simulation 
of  Z by  simply  obtaining  the  smallest  from  among  the 
realization  of  X^  X2,  and  X3-  The  advantage  of  the  pre- 
vious method,  based  on  (2.6)  for  instance,  is  that  only  one 
realization — that  of  an  exponential  in  that  specific  example 
— leads  to  the  realization  of  Z.  This  is  not  only  computa- 
tionally attractive,  but  seems  to  facilitate  the  application 
of  such  Monte  Carlo  variance  reduction  techniques  as  control 
and  antithetic  variables,  cf.,  Hammersley  and  Handscomb  [5]. 
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III.  OBTAINING  SPECIFIED  HAZARD  BEHAVIOR 


BY  SIMPLE  SAMPLING 


The  development  of  the  last  section  illustrates  one 
manner  in  which  hazard  behavior  may  be  conveniently  repre- 
sented and  simulated.  We  now  show  how  such  behavior  may 
alternatively  be  obtained  by  simple  simulation,  i.e.  from 
one  realization  of  a basic  (possibly  exponential)  random 
variable . 

Refer  to  (2.5) , in  which 


Z - G (X)  (3.1) 

and,  if  G(*)  is  monotonically  increasing. 


z(p)  » G(x(p) ) , # (3.2) 

z(p)  and  x(p)  being  the  p*100%  percentiles  of  Z and  X, 
respectively.  Then  the  counterpart  to  (2.23)  that 'results 
from  differentiation  of  (3.2)  is  the  expression 


Vz(P>>  - hx(x(p)>  G.,^p)  )■  > hxU(p))  (3^7^-  (3.3) 

Consequently,  if  one  specifies  h (z)  as  a suitable  func- 

z 


tion  of  the  "time"  z,  and  specifies  the  distribution  of 
the  stochastic  variable  X — and  hence  its  hazard,  hx — th€ 
results  a differential  equation  for  z(x)  = G(x) : 


h2(z) 


dz 

33? 


hxU) 


(3.4) 


integration  then  provides  the  desired  transformation,  G.  In 
other  words,  we  seek  z(x)  satisfying 
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(3.5) 


z x 

/ h (u)du  - / h (v)dv  , 


which  can  sometimes  be  carried  out  in  a useful  closed  form. 

Example  3.1.  Refer  to  the  example  of  Section  II,  wherein 

h is  given  by  the  expression  (2.50)  and  we  assume  that  X 
z 

is  exponential,  so  hx  is  constant.  Then  in  order  to  deter- 
mine G(x)  = z(x),  solve  the  equation 


A 

u + cx 


+ Bu  + X du  * / dv 


A ln(l  + j)  + jz^  + \z»x 


(3.6) 


Closed-form  solution  of  this  expression  for  z in  terms  of 
x is  of  course  impossible.  One  possible  approach  is  purely 
numerical:  find  an  approximate  solution,  zQ(x),  e.g.  the 

appropriate  solution  of  the  quadratic 


| z^  + \z  ■ X 


(3.7) 


and  then  correct  the  result  by  a few  Newton- Raphson  itera- 
tions. In  other  words,  put 


z1  (x) 


- X + V\  + 2Bx 


(3.8) 


now  apply  Newton  to  obtain  an  improved  solution 


A in  (1  + z^/a) 

z2  (x)  * zi  A77i7^rr+Tzpr 


(3.9) 


A 


The  process  can  be 


! 

! 


I 


which  will  be  feasible  if  0 < 
iterated  (the  numerator  will  change  after  the  first  iter- 
ation) . If  one  wishes  to  use  this  model  it  may  actually  be 
desirable  to  start  by  solving 


| z2  + (X  + £)z  - x * 0 
2 a 


(3.10) 


for  z^,  in  which  case  the  numerator  will  not  be  as  shown  in 
(3.9) ; convergence  may  be  more  rapid. 


Example  3.2.  Change  the  hazard  representation  of  the  previous 
example  as  follows:  let 


h„  (u) 


(v  + a) 


+ Bu  + X;  (A,  a,  B,  X > 0)  (3.11) 


then 


/ V">du  - a)  * I *2  + <3 


.12) 


Now  it  is  necessary  to  solve 


a ( z +V  + ! z2  + Az  = 55  ' 


(3.13) 


i.e. , the  cubic 


| z3  fa  | + X]  z2  + [|  + aX  - x]z  - ax  « 0 , (3.14) 

which  can  be  carried  out,  at  least  formally,  in  closed  form. 
Cnee  again  an  iterative  solution  that  begins  by  dropping  the 
cubic  term,  solving  the  resulting  quadratic  for  z^x),  and 
then  continuing  along  the  Newton- Raphson  road  may  be> 
successful.  Further  investigations  of  these  ideas  should 
be  conducted. 


IV.  COMPUTER  SIMULATION  AND  ESTIMATION  PROCEDURE 


A.  SIMULATION  AND  NUMERICAL  RESULTS 

In  previous  chapters  an  analytical  model  was  described 
for  the  failure  rate  function;  useful  formulas  were  also 
derived  from  the  model  (2.39)  and  (2.40).  Before  the  model 
is  used  in  realistic  situations,  it  will  be  convenient  to 
build  a computer  simulation  model  for  model  validation. 

1.  Algorithm 

First  a very  basic  simulation  model  was  built  for 
determining  the  general  shape  of  the  failure  rate  function 
associated  with  parameters  a,  8 and  X.  In  the  simulation 
model,  u was  selected  to  be  1.0  and  2.0,  8 was  selected  to 
be  0.05,  0.01,  0.005,  0.001  and  X was  selected  to  be  0.1. 
These  values  were  picked  arbitrarily,  it  is  only  stipulated 
that  a is  always  greater  than  8*  The  general  algorithm 
of  the  simulation  model  is  shown  in  Fig.  2. 

The  hazard  function  is  calculated  according  to  the 
model  (2.37)  and  the  system  logic  function  (2.40).  The 
results  were  shown  in  Fig.  3 and  Fig.  4. 

2 . Some  Comments 

These  simulation  results  show  that: 
a.  Parameter  a is  effective  when  z values  relatively 
have  small  values.  That  is,  it  influences  the  early 
failure  period. 
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INITIAL  VALUES  OF 
PARAMETERS 


GENERATING  RANDOM  NUM. 

FROM  EXPONENTIAL 

DISTRIBUTION 

GETTING  Z 

Z=G(X) 

VALUES 

GETTING  FAILURE  RATE 

FUNCTION 



PLOTTING  Z vs  FAILURE 
RATE  FUNCTION 


Fig.  2.  Basic  Algorithm. 
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Fig.  3.  The  Shape  or  Hazard  Function  with  Different  Values 


b.  Parameter  S is  effective  on  the  relatively  bigger  z 
values.  That  is,  it  describes  wearout  failures. 

c.  Parameter  X has  little  effect  on  the  shape  of  the 
curve;  it  is  a scale  factor. 

d.  The  last  important  observation  is  that  the  z values 
are  limited  by  the  parameter  0 such  that: 

2 - J 

When  z equals  1/8,  h (z)  goes  to  infinity. 

z 

B.  EMPIRICAL  DENSITY  FUNCTION  AND  ESTIMATION  PROCEDURE 

In  this  section,  an  estimation  procedure  for  parameters 
a,  8,  X is  defined  and  the  procedure  of  finding  the 
empirical  density  function  is  described. 

1.  Empirical  Density  Function 

Suppose  that  N^,  N,  A and  Nq  are  defined  as 


follows 


* the  frequency  of  data  points  for  each  time  interval 
between  a^  and  b^  for  i equals  1,2,3, ...,k. 

N » the  total  number  of  failures  where, 


N - l N. 
i-1  1 

the  length  of  each  interval 


A = b . - a . 
i i 
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NQ  ■ the  total  number  of  survivors  by  the  beginning  of 
each  interval  which  can  be  shown  as  follows: 


i-1 

N - I N 

j-1  3 


if  i * 1 


if  i = 2 , 3 , . . . ,k 


where  £ ",  N.  is  the  number  of  failures  before 
3“  3 i»l 

interval  i,  or  N - £ ^ N ^ is  the  number  of 


failures  after  interval  i. 


If  h2(z)  is  defined  as  the  density  function  of 
hazard,  then  the  relationship  between  the  frequences  of  the 
failure  data  and  density  function  of  hazard  can  be  given 
as  follows: 


£ Nq  • / hz(z)dz 


(4.1) 


where  i * 1,2 , 3, . . . ,k  and  a^  * (i-1) A,  b^  » id.  At  this 
point  some  approximation  can  be  made  in  expression  (4.1), 
such  that 


/ hz(z)dz  ~ (bi  - ai)*hz(i(bi  - ai) ) 


(4.2) 


for  small  values  of  A * b-a.  Then  this  approximation  (4.2) 
is  substituted  in  the  expression  (4.1),  it  turns  out  as 
follows : 


N.  * A*Nrt  hfiA) 
0 z 


(4.3) 
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« 


^z(lA)  " nTTI  ' 


X ~ 1 / 2 r « • • 


(4.4) 


for  the  general  case,  it  will  be; 


N. 

i 

N 


for  i = 1 


hz(iA) 


(4.5) 


N, 


i-1 


(N-  l N.) 

j-1  3 


for 


2,3, 


The  expression  (4.5)  will  be  used  in  the  calculation  of 
an  ampirical  hazard  function  and  also  used  in  the  proposed 
estimation  procedure. 

2 . Estimation  Procedure  for  Model  Parameters 

Two  approaches  can  be  used  for  this  problem.  The 
first  idea  is  to  use  the  relationship  between  a sample  pth 
percentile  and  the  related  probability  of  the  pth  percentile. 
A subsequent  idea  is  to  approach  the  problem  as  a nonlinear 
least  square  estimation  problem  for  a,  8,  X:  pick  a,  8 
X so  that  the  values  obtained  minimize  the  sum  of  squared 
errors  in  an  objective  function. 

a.  3-Percentile  Approach 

Suppose  Fx(x)  is  the  cumulative  probability 
function  of  the  exponential  distribution.  The  pth  percentile 
x(p)  is  equal  to  the  value  such  that: 

p - Fx(x(p))  * 1 - e"Xx(p)  (4.6) 
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Then 


x(p)  = F ^(p)  = - y £n(l-p)  (4.7) 

whe  re  X > 0 . 

If  z(p)  is  defined  to  be  the  pth  percentile 
in  the  bathtub  model  (2.37)  then  z(p)  and  x(p)  have  a 
definite  relationship  with  expression  such  that: 

2,  . 

, v ax  (p) 

2(P>  = (1  + axfprr  (1  + 3x(p)) 


by  substituting  (4.7) 


a{ - j In ( 1-p) 


z (p)  = 


jl  + a [-  j £n{l-p)  j | Jl  + 3 |"-  x In  (1-p) 


(4.8) 


or 


z(p) 


(a/X2) £2 (p) 

i1  + X 1 + f e(P)} 


(4.9) 


where  e(p)  » - In (1-p) . 

In  the  last  equation  (4.9)  , there  are  three  un- 
known variables  which  are  a,  3,  X.  If  three  independent 
equations  associated  with  expression  (4.9)  are  defined, 
clearly,  they  can  be  solved  for  the  unknown  parameters. 
Actually,  this  can  be  easily  done,  for  using  different  values 
of  percentiles  such  that: 
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z(Px) 


(a/A  ) 


I1  + y£  (piJ}  { 1 + I £(pi)} 


z(P2) 


z(P3) 


(a/A2)e2(P2) 


i1  +f  e(P2}}  I1  + X £(P2>F 

(a/A2)  e2  (P3) 

{1*fS(P3,}{’-+I£<P3>} 


where  z(P^)  and  e (P^)  are  the  known  values  associated 
with  P^.  Actually,  to  get  the  percentile  values,  two  kinds 
of  approaches  can  be  made.  First  they  can  be  computed 
directly  from  data  using  the  simple  statistics  method,  such 
as  z(jj  is  (approximately)  z(i/(n+l)).  Second  they  can 
be  computed  from  the  empirical  density  function.  Generally 
the  choice  depends  on  the  form  of  the  data  that  are  avail- 
able. 

To  decide  for  the  effectiveness  of  this  type  of 
estimation,  another  basic  computer  simulation  is  made  by 
modifying  the  first  computer  simulation  algorithm.  This 
algorithm  is  shown  in  Fig.  5. 

The  simulation  was  run  four  hundred  times  and 
in  each  run  the  sample  size  was  assumed  to  be  n = 50. 
Initially  the  parameters  a,  B and  A were  taken  to  be, 
respectively,  1.0,  0.05  and  0.1.  The  percentiles  P^,  P2 
and  P^  were  used  in  each  run  such  that  0.1,  0.5  and  0.9 
respectively.  The  estimation  results  are  shown  in  Table  I. 
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SELECTION  OF  DIFFRENT 
VALUES  OF  PARAMETERS 


GENERATING  RANDOM  NUM 
FROM  EXPONENTIAL 
DISTRIBUTION 


GET  THE  FREQUENCY 
DISTRIBUTION  FOR 
CERTAIN  INTERVAL 


FIND  EMPIRICAL  DENSITY 
FOR  HAZARD, ESTIMATE 
THE  PAREMETERS 


Fig . 5 .Simulation  algoritm  for  estimation  procedure 
and  emprical  density  function. 
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As  can  be  seen,  the  estimated  parameters  a,  0, 

X have  quite  different  values  in  all  runs.  To  judge  their 
sampling  distributions,  histograms  were  drawn  separately; 
Figures  6 , 7 and  8 are  the  histograms  of  the  estimated 
values  of  a,  8,  X,  respectively. 

It  is  suggested  by  Figures  6 and  7 for  X and  0, 
that  normality  may  be  assumed  with  sample  means  0.9246, 
0.05488  and  sample  variances  0.001065,  0.0002246,  respec- 
tively.  Their  mean,  median,  trimean  and  midmean  are  pretty 
close  to  each  other.  But,  in  the  case  of  a,  the  histogram 
(Fig.  8)  is  a quite  different  picture  which  looks  something 
like  an  exponential  distribution  instead  of  normal  distri- 
bution. The  reason  is  clear,  because  the  negative  values 
of  estimated  a were  not  taken  into  account.  They  are 
physically  infeasible. 

b.  Nonlinear  Least  Squares  Approach 

The  least  squares  criterion  used  here  can  be 
stated  formally  and  generally  as  follows; 

N - 2 

minimize  £ (Y.  - Y.)  (4.10) 

i-1  1 1 

where  N is  the  number  of  observations,  Y^  is  the  fitted 
value  of  Y^.  In  this  case,  expression  (4.10)  can  be  re- 
written as  follows 

N - 

min  l {Z.  - Z , (a , 0 , X)  } (4.11) 

i-1  1 1 
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TABLE  I 


COMPUTER  SIMULATION  RESULTS  FOR  ESTIMATION  PROCEDURE 
Sample  size  n = 50;  Number  of  Total  Runs  m=  400; 


Value  of  parameters  a,  0,  X Introduced:  a-1,  0-0.05,  X-0.1. 


Number  of 
Runs 

X 

a 

0 

1 

0.0858 

2.2255 

0.0549 

2 

0.06591 

1.2745 

0.0779 

3 

0.08381 

9.498 

0.0688 

4 

0.0881 

0.9304 

0.0  39 

5 

0.1286 

3.7429 

0.0498 

6 

0.0883 

2.114 

0.069 

7 

0.0785 

2.8281 

0.0593 

8 

0.0945 

0.4073 

0.0584 

9 

0.089 

0.3833 

0.0452 

10 

0.1344 

1.3351 

0.0426 

11 

0.0677 

0.2111 

0.0729 

12 

0.1304 

2.4871 

0.0512 

13 

0.0714 

0.3211 

0.0623 

14 

0.069 

0.3547 

0.0657 

15 

0.0775 

0.4813 

0.0772 

16 

0.1618 

2.7416 

0.027 

17 

0.0827 

0.927 

0.0705 

18 

0.0816 

1.1663 

0.0836 

19 

0.0515 

0.3777 

0.0775 

20 

0.1089 

1.6439 

0.0522 

21 

0.0767 

0.8667 

0.0605 

22 

0.0499 

0.3414 

0.0748 

23 

0.0569 

0.2118 

0.0608 

24 

0.0648 

0.2689 

0.0715 

25 

0.1415 

2.1178 

0.0518 

26 

0.1636 

1.3635 

0.0186 

27 

0.0947 

6.101 

0.0664 

TABLE  I Cont 


Number  of 
Runs 

\ 

a 

0 

28 

0.1063 

0.5235 

0.0623 

29 

0.0825 

2.6602 

0.059 

30 

0.13 

0.3412 

0.0303 

31 

0.817 

0.24 

0.0685 

32 

0.0843 

0.4377 

0.0393 

33 

0.1445 

4.6941 

0.0442 

34 

0.0737 

0.2956 

0.0665 

35 

0.045 

0.1372 

0.0646 

36 

0.1108 

1.1026 

0.0427 

37 

0.1277 

1.091 

0.052 

38 

0.0514 

0.1362 

0.062 

39 

0.1165 

0.6  79  4 

0.0469 

40 

0.0954 

3.7845 

0.0587 

41 

0.0563 

0 .4123 

0.0738 

42 

0.0684 

0.2903 

0.0574 

43 

0.1263 

2.059 

0.0357 

44 

0.1335 

7. 1595 

0.0226 

45 

0.0708 

0.465 

0.063 

46 

0.0872 

2. 164 

0.0657 

47 

0.0834 

0.3568 

0.0567 

48 

0.0429 

0.127 

0.0694 

49 

0.0428 

0.1295 

0.0639 

50 

0.0834 

0.6807 

0.0489 

51 

0.1141 

1.5349 

0.0566 

52 

0.1057 

1.1757 

0.0441 

53 

0.1104 

5.0692 

0.0433 

54 

0.1426 

0.8845 

0.0244 

55 

0.10  84 

0.5775 

0.0505 

56 

0.1048 

1. 8399 

0.0562 

57 

0.0561 

0.3197 

0.0643 

58 

0.0607 

0.4564 

0.0607 

59 

0.0724 

0.3133 

0.0551 

43 


TABLE  I Cont 


Number  of 
Runs 

A 

a 

0 

60 

0.0957 

0.9584 

0.0557 

61 

0.0804 

0.7378 

0.0423 

62 

0.0888 

0.1718 

0.0393 

63 

0.064 

0.3779 

0.9543 

64 

0.0324 

0.2202 

0.0791 

65 

0.0695 

0.2428 

0.0564 

66 

0.0903 

6.4194 

0.0542 

67 

0.0741 

0.4118 

0.0718 

68 

0.0462 

0.167 

0.0671 

69 

0.0805 

1.0279 

0.0662 

70 

0 .0  469 

0.2157 

0.0693 

71 

0.0487 

0.2457 

0.0688 

72 

0.0682 

1.1847 

0.0665 

73 

0.0538 

0.1964 

0.0573 

74 

0.1516 

2.099 

0.0352 

75 

0.0  76  3 

0.3067 

0.0679 

76 

0.0558 

1.5589 

0.0  79 

77 

0.0463 

0.2252 

0.0612 

78 

0.095 

5.6243 

0.0627 

79 

0.190  3 

6.0 

0.0604 

80 

0.0822 

0.9989 

0.0738 

81 

0.066 

0.3422 

0.0926 

82 

0.0784 

0.3336 

0.0517 

83 

0.0499 

0.1453 

0.062 

84 

0.0667 

0.1357 

0.064 

85 

0.0431 

0.0761 

0.0568 

86 

0.1485 

6.1674 

0.0262 

87 

0.1154 

1.0451 

0.0476 

88 

0.0043 

0.349  3 

0 0648 

89 

0.0983 

1.3101 

0.0422 

90 

0.138 

6.7408 

0.0443 

TABLE  I Cont 


Number  of 
Runs 

X 

a 

6 

91 

0.1117 

0.7281 

0.058 

92 

0.065 

0.4284 

0.0482 

93 

0.0638 

1.2754 

0.0611 

94 

0.1603 

1.003 

0.0302 

95 

0.0436 

0.1647 

0.085 

96 

0.0905 

6.2455 

0.0703 

97 

0.1309 

1.4348 

0.0431 

98 

0.0787 

0.5245 

0.0686 

99 

0.1927 

7.801 

0.009 

100 

0.0683 

0.525 

0.0698 

101 

0.1163 

1.7146 

0.0355 

102 

0.1005 

1.5795 

0.0463 

103 

0.1492 

1.5169 

0.0465 

10  4 

0.1511 

2.7243 

0.318 

105 

0.0823 

0.16  34 

0.0438 

106 

0.1493 

6.9467 

0.0437 

107 

0.1185 

9 . 4492 

0.0603 

108 

0.1039 

0.589 

0.0406 

109 

0.0832 

0.39  35 

0.0612 

no 

0.059 

0.145 

0.0508 

111 

0.133 

8.6071 

0.0  39  7 

112 

0.0877 

0.4751 

0.0448 

113 

0.0772 

0.2752 

0.0624 

114 

0.0659 

0.4782 

0.0583 

115 

0.0796 

0.2978 

0.0537 

116 

0.1128 

1.668 

0.0601 

117 

0.0974 

1.2931 

0.0515 

118 

0.1105 

0.3394 

0.0468 

119 

0.1216 

1.0529 

0.0655 

120 

0.1046 

2.0384 

0.0591 

121 

0.128 

3.2468 

0.0405 

TABLE  I Cont. 


Number  of 
Runs 

X 

01 

8 

122 

0.0818 

1.092 

0.0584 

12  3 

0.0924 

2.4295 

0.0439 

124 

0.0732 

0.1916 

0.0495 

125 

0.1284 

1.8515 

0.0407 

126 

0.1376 

2.6651 

0.0282 

12  7 

0.0877 

7.8429 

0.0703 

128 

0.0865 

1.1942 

0.0438 

129 

0.1404 

2.8991 

0.0326 

130 

0.0452 

0.0978 

0.0663 

131 

0.0928 

3.0463 

0.0471 

132 

0.0668 

0.4602 

0.0551 

133 

0.2089 

0.7373 

0.0119 

134 

0.0946 

0.4011 

0.038 

135 

0.1352 

1.2162 

0.0284 

136 

0.0597 

0.2241 

0.0539 

137 

0.106 

1.7821 

0.0354 

138 

0.0508 

0.2689 

0.0596 

139 

0.0522 

0.1244 

0.055 

140 

0.104 

0.8464 

0.0603 

141 

0.1002 

0.817 

0.0711 

142 

0.0692 

1.4214 

0.0744 

143 

0.0995 

0.5442 

0.0435 

144 

0.0697 

0.4966 

0.0487 

145 

0.0737 

0.8995 

0.0708 

146 

0.059 

0.2998 

0.0654 

147 

0.1235 

1.1251 

0.045 

148 

0.1023 

0.615 

0.0349 

149 

0.1 

0.9087 

0.0466 

150 

0.1242 

2.9267 

0.0356 

151 

0.0966 

0.4989 

0.0616 

152 

0.0913 

2.5416 

0.0638 

153 

0.1301 

1. 3897 

0.0434 

TABLE  I Cont 


Number  of 
Runs 

\ 

a 

6 

154 

0.1033 

2.0535 

0.0466 

15‘> 

0.0638 

0.2219 

0.0594 

156 

0.1408 

3.6331 

0.0434 

157 

0.1311 

4.9526 

0.0495 

158 

0.0608 

0.469 

0.0462 

159 

0.0574 

0.5886 

0.0723 

160 

0.0998 

4.912 

0.0554 

161 

0.1114 

0.4898 

0.0363 

162 

0.0415 

0.1827 

0.0854 

16  3 

0.1586 

7.5257 

0.0  39  7 

164 

0,0723 

0.9055 

0.0637 

165 

0.0984 

2.1521 

0.0675 

166 

0.1011 

2.1535 

0.0573 

16  7 

0.093 

4.8488 

0.0555 

168 

0.1141 

4.2636 

0.0443 

169 

0.1953 

4.6512 

0.0452 

170 

0.0568 

0.8364 

0.07099 

171 

0.0927 

5.130  3 

0.068 

172 

0.0632 

0.7399 

0.069 

173 

0.0536 

0.2732 

0.07138 

174 

0.1033 

0.9515 

0.0787 

175 

0.0633 

1.6597 

0.0791 

176 

0.1414 

1.2547 

0.0318 

177 

0.1089 

4.1786 

0.051 

178 

0.1165 

0.5441 

0.0432 

179 

0.1731 

4.2913 

0.0192 

180 

0.0679 

0.1064 

0.0451 

181 

0.1682 

3.1417 

0.0332 

182 

0.0868 

0.6399 

0.0685 

183 

0.0807 

1.046 

0.0607 

184 

0.0778 

1.185 

0.0694 

47 


Number  of 
Runs 


at 


0 


185 

0.0769 

0.1999 

0.0009 

186 

0.0396 

0.1194 

0.064 

187 

0 . 10 12 

0.499 

0.485 

188 

0.1389 

9.70  85 

0.0334 

189 

0.0665 

1.1661 

0.0599 

190 

0.0471 

0.1374 

0.0643 

191 

0.1316 

4.1676 

0.0259 

192 

0.0687 

0.6692 

0.0928 

193 

0.0866 

0.5419 

0.0484 

194 

0.0966 

1.579 

0.0524 

195 

0.0773 

1.0102 

0.0668 

196 

0.1076 

0.7413 

0.0704 

197 

0.0857 

0.8715 

0.0548 

19  8 

0.0803 

1.110 

0.0549 

199 

0.052 

0.1396 

0.0598 

200 

0.0757 

0.9558 

0.0622 

cell 


FIG.  6.  Histogram  for 


counts 


The  goal  is  to  choose  values  a,  8 and  X which  minimize 
the  expression  (4.11).  In  the  linear  case,  this  can  be  done 
using  elementary  calculus  by  taking  the  partial  derivatives 
with  respect  to  a,  8 and  X , setting  each  of  them  equal 
to  zero  and  solving  the  resulting  linear  equations  (Normal 
equations) . But  the  present  problem  is  different  because 
the  expression  (4.11)  is  highly  nonlinear  in  the  parameters, 
and  indefinite  in  terms  of  convexity.  This  is  why  a non- 
linear programming  technique  is  used. 

There  are  a few  general  approaches  to  the  solu- 
tion of  the  nonlinear  estimation  problem.  One  of  them  is 
the  direct  optimization  approach. 


Specifically,  equation  (4.9)  is  considered  as 


follows : 


P1P2z  (P) 


+ P2e 


+ P3e 


(4.12) 


where 


In  equation  (4.12),  z(p)  and  e (p)  are  known  values, 
are  unknown  variables.  If  it  is  rewritten  as  the  sum  of 
squared  errors  the  objective  function  will  be: 


N | 

l 

i-1  ( 


- 


°lp2e  (Pi* 


+ P2e 


+ p3e(Pi))  j 


(4.13) 


Now  nonlinear  optimization  problems  can  be  stated  such  that 


N ( 

Min  l z (P. ) - 

i-M 


plp2e 


+ p3e 


Subject  to 


P1  - p2  - p3  - 0 


In  order  to  perform  the  optimization,  i.e.  to 
solve  this  nonlinear  optimization  problem,  the  GRG  package 
(Generalized  Reduced  Gradient)  [6]  was  used;  it  is  very 
convenient  for  this  case.  Ten  runs  were  made  with  five 
different  initial  points,  both  with  analytically  computed 
derivatives  and  with  numerically  computed  derivatives.  The 
first  four  initial  points  were  selected  arbitrarily.  The 
last  one  was  picked  such  that: 


(4.14) 


where  2(u  is  t^ie  smallest  number  in  the  sample,  2v5q) 
is  the  largest  number  in  the  sample,  and  X is  the  natural 
logarithm  of  p = 0.5  divided  by  the  median  of  the  sample. 
Initial  points  are  shown  in  Table  II. 

The  optimization  results  of  a,  3,  and  p 
associated  with  p^,  p2  and  p^  are  tabulated  in 
Table  III  and  Table  IV,  respectively. 
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As  it  is  seen  from  Tables  III  and  IV,  the  program 


i 

( 


! 


was  terminated  at  different  optimal  points.  Because  the 
objective  function  is  highly  nonlinear  and  apparently  in- 
definite (neither  convex  nor  concave) , it  may  sometimes 
have  stopped  at  a local  minimum  rather  than  at  a global 
optimum  point.  But  if  the  results  are  compared  to  the 
3-percentile  results,  the  nonlinear  least  square  estimates 
show  much  greater  accuracy  and  seem  more  consistent. 


TABLE  II 

INITIAL  GUESS  POINTS 


# of 
Points 

P1 

P2 

p3 

I 

50.0 

2.0 

30.0 

II 

10.0 

10.0 

0.5 

III 

80.0 

80.0 

5.0 

IV 

1.0 

1.0 

1.0 
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TABLE  IV 


Actual  values  X * 0.1,  a * 1.0,  8 = 0.05 


of 

ins 

X 

a 

8 

Point 

A 

B 

A 

B 

A 

B 

I 

0.1202 

0.1203 

343.512 

676.800 

0.0400 

0.0400 

II 

0.1199 

0.1203 

129.591 

968.150 

0.0401 

0.0401. 

1 

III 

0.1168 

0.1203 

8.898 

1506.926 

0.0409 

0.0400 

IV 

0.1158 

0.1204 

6.577 

874.957 

0.0412 

0.0385 

V 

0.1203 

0.1201 

4.182 

955.629 

0.0400 

0.0401 

I 

0.0984 

0.0863 

1 

0.0478 

0.0510 

II 

0.0865 

0.0863 

0.0509 

0.0510 

2 

III 

0.0977 

0.0  86  3 

6.158 

1.297 

0.0480 

0.0510 

IV 

0.0798 

0.0863 

0.810 

1.292 

0.0527 

0.0511 

V 

0.0500 

0.0  86  3 

0.154 

1.296 

0.0583 

0.0511 

I 

0.1390 

0. 1434 

8.123 

392.550 

0.0343 

0.0331 

II 

0.1430 

0.1434 

138.531 

1520.067 

0.0332 

0.0332 

3 

III 

0.1440 

0.1436 

204.078 

1465.586 

0.0328 

0.0330 

IV 

0.1434 

0.1433 

414.903 

1177.650 

0.0331 

0.0332 

V 

0.0438 

0.1434 

0.055 

2056.280 

0.0515 

0.0331 

I 

0.0878 

0.0843 

0.118 

0.10  7 

0.0394 

0.0397 

II 

0.0890 

0.0873 

0.12  3 

0.116 

0.0392 

0.0  39  3 

4 

III 

0.0891 

0.0856 

0. 12  3 

0.111 

0.0393 

0.0395 

IV 

0.0870 

0.0829 

0.115 

0.10  3 

0.0  39  4 

0.0  39  8 

V 

0.0899 

0.0860 

0.125 

0.112 

0.0  390 

0.0394 

I 

0.0507 

0.0508 

4.944 

0.509 

0.0530 

0.0775 

II 

0.0695 

0.0509 

0.847 

0.513 

0.0988 

0.0775 

5 

III 

0.0692 

0.0508 

3.882 

0.508 

0.0749 

0.0775 

IV 

0.0540 

0.0508 

0.826 

0.509 

0.0717 

0.9775 

V 

0.0676 

0.0507 

0.6  72 

0.506 

0.1034 

0.0776 
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TABLE  IV  Cont 


# of 
Runs 

Initial 

Point 

X 

a 

8 

A 

B 

A 

B 

A 

B 

I 

0.1384 

0.0994 

9.072 

■ 

0.0494 

0.0591 

II 

0.1019 

0.0996 

0.628 

0.0585 

0.0590 

6 

III 

0.1363 

0.0998 

6.105 

0.0499 

0.0589 

IV 

0 . 10  79 

0.0998 

0.810 

0.582 

0.0571 

0.0590 

V 

0.1272 

0.0994 

2.439 

0.571 

0.0525 

0.0590 

I 

0.1158 

0.0738 

6.780 

0.301 

0.0405 

0.0502 

II 

0.0755 

0.0738 

0.324 

0.301 

0.0496 

0.0502 

7 

III 

0.0752 

0.0743 

0.319 

0.309 

0.0499 

0.0501 

IV 

0.0757 

0.0739 

0.328 

0. 30  8 

0.0499 

0.0501 

V 

0.1188 

0.0738 

28.839 

0.301 

0.0397 

0.0502 

I 

0.0753 

0.0748 

7.325 

5.783 

0.0645 

0.0631 

II 

0.0739 

0.0740 

4.647 

4.763 

0.0634 

0.0634 

8 

III 

0.0735 

0.0735 

4.219 

4.221 

0.0635 

0.0635 

IV 

0.0737 

0.0741 

4.361 

4.831 

0.0634 

0.0633 

V 

0.0170 

0.0010 

0.033 

0.0724 

0.0634 

I 

0.1257 

0.1252 

11.646 

9.302 

0.0275 

0.0276 

II 

0.1252 

0.125  3 

10.181 

9.374 

0.0275 

0.0275 

9 

III 

0.1240 

0.1246 

6.983 

8.598 

0.0278 

0.0278 

IV 

0.1254 

0.1254 

11.010 

10.296 

0.0276 

0.0275 

V 

0.1017 

0.1255 

0.754 

10.407 

0.0330 

0.0275 

I 

0.1387 

0.1438 

5.688 

25.170 

0.0326 

0.0313 

II 

0.1441 

0.1428 

28.830 

29.290 

0.0312 

0.0310 

10 

III 

0.1400 

0.1440 

6.290 

27.747 

0.0326 

0.0313 

IV 

0.1437 

0.1441 

21.370 

26.590 

0.0313 

0.0312 

V 

0.1415 

0.1427 

10.0  30 

14.610 

0.0310 

0.0316 

C.  TEST  PROCEDURE  FOR  PARAMETERS 

In  order  to  determine  whether  the  model  reasonably  fits 
the  data  a simple  test  procedure  is  applied.  It  is  to 
exhibit  the  simple  comparison  of  the  actual  points  z(P) 
and  the  estimated  values  z(P),  using  the  estimated 
parameters. 

Basically  z(P)  can  be  obtained  as  follows: 


£ Hn(l-p) 

A 


(4.15) 


Then  if  it  is  substituted  in  expression  (2.37),  z (p)  will 


A AO 

axz(P) 

(1  + ax (P) ) (1  + 0x(P) ) 


(4.16) 


In  equation  (4.16),  all  variables  are  known  so  estimated 
quantiles  can  be  easily  obtained. 

In  Fig.  9 actual  values  and  estimated  values  are  plotted. 
Also  Fig.  10  displays  residuals,  which  are  the  difference 
between  actual  and  estimated  values. 
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V.  NUMERICAL  APPLICATIONS  TO  TWO  SETS 


OF  REAL  DATA 


In  this  section  the  two  failure  data  sets  studied  come 
from  these  sources:  Oral  irrigator  [7]  and  human  life  [8]. 
They  are  used  for  numerical  examples  in  that  these  data  are 
fitted  using  the  model  (2.6). 


A.  ORAL  IRRIGATOR 

The  data  used  was  obtained  from  the  Commun.  Statist. - 
Theor.  meth.,  Colvert  and  Bordman  [7],  The  data  was  col- 
lected such  that  100  oral  irrigators  were  placed  on  the  test 
was  terminated  at  700  time  units.  During  the  life  test,  98 
oral  irrigators  failed;  another  2 oral  irrigators  survived. 
The  ordered  observed  times  to  failure  of  the  oral  irrigators 
is  tabulated  in  Table  V. 

Using  the  data  of  Table  V,  the  parameters  a,  8,  and  X 
are  estimated  by  means  of  the  3-percentile  approach.  To 
demonstrate  the  differences  between  estimated  values, 
various  different  combinations  of  p values  w*jre  used.  The 
results  are  shown  in  Table  VI. 

Examination  of  Table  VI  indicates  that  the  3-percentile 
approach  may  produce  estimates  having  great  differences  for 
different  percentile  values,  except  here  in  the  case  of  the 
first  two  combinations.  Also  the  first  two  estimations  seem 
to  have  accepatable  limiting  age.  Moreoever,  the  last  two 
combinations  do  not  fit  well,  as  judged  from  the  residuals. 


TABLE  V 


TIME  TO  FAILURE  OF  ORAL  IRRIGATORS 


1.75 

7.02 

7.58 

9.76 

15.02 

15.57 

17.39 

19.55 

22.47 

23.24 

23.96 

25.05 

32.44 

36.87 

42.76 

43.14 

46.95 

56.33 

58.99 

59.08 

60.37 

61.01 

77.86 

86.45 

88.50 

103.06 

104.34 

105.85 

117.46 

120.11 

122.28 

122.61 

129.31 

130.42 

137.57 

142.27 

142.98 

148.29 

150.79 

151.21 

155.62 

157.93 

160.72 

169.79 

186.26 

197.60 

224.83 

233.64 

242.07 

256.86 

260.77 

261.68 

277.99 

283.95 

288.94 

295.48 

314.76 

316.06 

332.07 

339.46 

362.61 

369.47 

370.74 

491.06 

40  3.39 

414.78 

426.71 

459.62 

455.84 

457.94 

466.61 

468.64 

469.09 

476.42 

481.41 

481.82 

488. 15 

490.06 

493.67 

494.38 

503.72 

508.93 

509.01 

418.32 

532.29 

534.62 

545.23 

547.41 

558. 41 

571.10 

585.52 

539.11 

592.93 

60  7.15 

623.15 

647.91 

TABLE  VI 


ESTIMATED  VALUES  OF  a , 3 , X BY  USING  3-PERCENTILE 


a 

P2 

P3 

A 

X 

a 

8 

1/6 

.i 

.5 

.9 

0.00153 

0.00634 

0.001028 

972.76 

in 

CM 

• 

.5 

.9 

0.00158 

0.00710 

0.001018 

982.31 

.1 

.5 

.75 

0.00239 

0.01835 

0.000228 

4385.96 

.25 

.5 

.75 

0.00273 

0.04514 

0.000077 

12987.01 

Next  the  nonlinear  optimization  approach  was  applied  to 
estimate  parameters  a,  $ and  X associated  with  p^,  p2 
and  p3,  both  using  the  derivative  and  without  derivative  by 
taking  into  account  9 8 ordered  values.  The  results  are 
tabulated  in  Tables  VII  and  VIII  for  p1#  p2,  P3  and  X 
a,  8 respectively.  Some  initial  points  in  Table  II  are 
used.  Table  VII  indicates  that  the  estimates  obtained  from 
the  nonlinear  estimation  method  are  more  consistent  than 
those  from  the  3-percentile  approach.  Also,  the  boundary 
points,  1/S,  seem  reasonable.  Finally,  comparison  of  the 
actual  values  of  z(p)  and  estimated  values  of  z(p)  indi- 
cate that  the  fitting  is  reasonable.  The  residuals  of 
fitted  model  and  the  comparison  of  the  actual  and  estimated 
values  are  plotted  in  Fig.  11  and  Fig.  12. 

B.  HUMAN  LIFE  (MORTALITY)  DATA 

The  data  used  in  this  example  was  obtained  from  the 
1969-71  life  table  [8]  for  white  females  in  the  United  States. 
The  table  has  been  prepared  from  a history  of  100,000  persons: 
the  number  of  surviving,  and  the  number  dying  has  been  given 
for  each  age  interval.  Look  at  Table  IX. 

1.  Estimation  by  Using  3-Percentile  Approach 

The  life  data  is  first  fitted  to  the  model  (2.37) 
by  using  the  3-percentile  approach.  Certain  pth  percentile 
values  are  selected  and  used  in  the  estimation  process. 

These  results  are  tabulated  in  Table  X.  Investigation  of 
the  last  eight  combinations  of  percentiles  in  these  tables 
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ESTIMATED  PARAMETERS  FOR  ORAL  IRRIGATORS 


indicates  that  the  limiting  age,  1/(3,  is  approximately  110, 
which  is  reasonable.  However,  the  first  combination  of 
percentiles  gives  a good  fit  in  terms  of  the  sum  of  the 
squared  errors.  But  for  these  limiting  age  is  reduced  to 
about  95  years.  This  apparently  means  that  the  fit  of  the 
model  does  not  represent  the  age  above  95.  Thus  it  can 
be  said  that  there  is  a trade-off  between  1/3  and  the  sum 
of  squared  errors.  The  present  model  simply  does  not  seem 
to  fit  the  mortality  data  very  well. 

2 . Estimation  By  Using  Nonlinear  Least  Squares,  and 
Application  of  a Nonlinear  Programming  Algorithm 
In  the  previous  section  we  discussed  the  fact  that 
one  of  the  problems  encountered  in  the  3-percentile  estima- 
tion procedure  was  the  high  uncertainty  of  fitting,  as 
measure  by  the  sum  of  squared  errors.  In  order  to  reduce 
this  variability,  a nonlinear  estimation  process  is  applied. 
Previously  it  was  noted  that  the  GRG  package  can  be  used 
with  analytically  computed  derivatives  and  without 
derivatives.  If  it  is  used  with  analytically  computed 
derivatives,  it  is  necessary  to  provide  another  subprogram 
by  the  user.  Otherwise  the  GRG  package  will  provide  the 
derivative,  computed  numerically  and  automatically. 

The  derivative  of  the  objective  function  (4.13)  is 
the  normal  equations  such  that 


TABLE  IX 

MORTALITY  DATA  FOR  WHITE  FEMALES,  1969-71 


Time 

Interval 

# of 
Survivors 

# of 
Dying 

Time 

Interval 

# of 

Survivors 

# of 
Dying 

0-  1 

100,000 

1532 

27-28 

97,165 

70 

1-  2 

99,468 

100 

28-29 

97,095 

73 

2-  3 

98,368 

65 

29-30 

97,022 

77 

3-  4 

98,303 

54 

30-31 

96,945 

81 

4-  5 

98,249 

46 

31-32 

96,864 

87 

5-  6 

98,203 

39 

32-33 

96,777 

93 

6-  7 

98,164 

35 

33-34 

96,684 

101 

7-  8 

98,129 

32 

34-35 

96,583 

109 

8-  9 

98,097 

29 

35-36 

96,474 

118 

9-10 

98,068 

26 

36-37 

96,356 

12  8 

10-11 

98,042 

24 

37-38 

96,228 

141 

11-12 

98,018 

24 

38-39 

96,087 

153 

12-13 

97,994 

25 

39-40 

95,932 

170 

13-14 

97,969 

30 

40-41 

95,762 

185 

14-15 

97,939 

37 

41-42 

95,577 

201 

15-16 

97,902 

45 

42-43 

95,376 

220 

16-17 

97,857 

54 

43-44 

95,156 

242 

17-18 

97,803 

60 

44-45 

94,914 

263 

18-19 

97,743 

62 

45-46 

94,649 

291 

19-20 

97,681 

63 

46-47 

94,358 

317 

20-21 

97,618 

63 

47-48 

94,041 

344 

21-22 

97,555 

63 

48-49 

93,697 

372 

22-23 

97,492 

63 

49-50 

93,325 

401 

23-24 

97,429 

64 

50-51 

92,924 

433 

24-25 

97,365 

66 

51-52 

92,491 

469 

25-26 

97,299 

66 

52-53 

92,022 

506 

26-27 

97,233 

68 

53-54 

91,516 

546 

71 


Time 

interval 

# of 

Survivors 

» of 
Dying 

Time 

Interval 

# of 

Survivors 

« of 
Dying 

54-55 

90,970 

587 

82-  83 

41,215 

3,586 

55-56 

90,383 

632 

83-  84 

32,629 

3,589 

56-57 

89,751 

680 

84-  85 

34,040 

3,550 

57-5  8 

89,071 

730 

85-  86 

30,490 

3,495 

58-59 

88,341 

781 

86-  87 

26,995 

3,425 

59-60 

87,560 

834 

87-  88 

23,570 

3,286 

60-61 

86,726 

911 

88-  89 

20,284 

3,072 

61-62 

85,835 

953 

89-  90 

17,212 

2,806 

62-63 

84,882 

1,022 

90-  91 

14,406 

2,531 

6 3-64 

83,860 

1,098 

91-  92 

11,875 

2,262 

64-65 

82,762 

1,183 

92-  93 

9,613 

1,982 

65-66 

81,579 

1,275 

93-  94 

7,631 

1,694 

66-67 

80,304 

1,375 

94-  95 

5,937 

1,411 

67-68 

78,929 

1,486 

95-  96 

4,526 

1,145 

68-69 

77,443 

1,607 

96-  97 

3,381 

905 

69-70 

75,836 

1,735 

97-  98 

2,476 

696 

70-71 

74,101 

1,862 

98-  99 

1,780 

524 

71-72 

72,239 

1,993 

99-100 

1,256 

384 

72-73 

70,246 

2,141 

100-101 

872 

277 

73-74 

68,105 

2,312 

101-102 

595 

195 

74-75 

65,793 

2,503 

102-103 

400 

135 

75-76 

63,290 

2,693 

10 3-104 

265 

92 

76-77 

60,597 

2,872 

104-105 

173 

61 

77-78 

57,725 

3,039 

105-106 

112 

41 

78-79 

54,686 

3,187 

106-107 

71 

26 

79-80 

51,499 

3,317 

107-108 

45 

17 

80-81 

48,182 

3,434 

108-109 

28 

11 

81-82 

44,748 

3,533 

109-110 

17 

6 

2 
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TABLE  X 


ESTIMATED  VALUES  OF  a,  0,  A FOR  MORTALITY  DATA  USING 

3-PERCENTILE  APPROACH 


Percentile 

Values 

22 

56 

92 

20 

56 

92 

17 

56 

92 

16 

56 

92 

17 

56 

109 

18 

56 

109 

19 

56 

109 

17 

56 

110 

18 

56 

110 

19 

56 

110 

20 

56 

110 

0.000674 

0.000550 

0.000508 

0.000440 

0.000862 

0.000886 

0.000908 

0.000876 

0.000899 

0.000920 

0.000941 


0.1190 

0.0404 

0.0302 

0.0199 

0.1676 

0.2389 

0.3708 

0.1805 

9.2613 

0.4188 

0.9070 


A 

8 

1 

A 

1/S 

0.010550 

94.7 

0.010520 

95.0 

0.010550 

94.7 

0.010550 

94.7 

0.009069 

110.2 

0.009068 

110.2 

0.009067 

110.2 

0.008991 

111.2 

0.008989 

111.2 

0.008988 

111.2 

0.008988 

111.2 

3 


r 
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N is  110  in  the  case  of  human  mortality  data. 

The  results  of  estimation  were  shown  in  Tables  XI 
and  XII  for  pi  and  a,  8,  A respectively.  Also  Fig.  13 
and  Fig.  14  demonstrate  the  actual  and  estimated  values  and 
residuals  of  fitting. 

Examination  of  Table  XI  and  Table  XII  indicates 
that  the  objective  function  values  (the  sum  of  the  squared 
errors)  of  the  fitting  are  smaller  than  the  objection  func- 
tion values  of  3-percentile  approach.  Also  the  estimated 


parameter  values  $ and  X associated  with  and  p ^ 

are  very  close  to  each  other.  But  estimated  values  of  a 
show  considerable  difference  for  the  various  initial  values. 
However,  there  are  two  important  facts  to  notice.  First, 
when  a increases  significantly,  the  objective  function 
values  remain  almost  constant.  That  is,  it  is  not  effected 
significantly  on  the  objective  function  values  in  this  case. 
Second,  this  analytical  model  (2.6)  does  not  represent 
deaths  beyond  the  age  95.  However,  modifications  can  be 
made  in  the  model  (2.6)  for  this  kind  of  difficulty.  Some 
ideas  will  be  discussed  in  a later  section. 

3.  Model  Modifications 

In  the  previous  section  it  was  noted  that  the  model 

(2.6)  has  relatively  great  errors  and  does  not  accurately 

describe  probability  of  death  at  age  greater  than  95  in  the 

human  life  example.  This  means  that  the  hazard  function 

h_(z)  increases  too  rapidly  in  the  wearout  period.  If  it 
z 

is  possible  to  slow  the  rate  of  increase  of  hazard  perhaps  a 
better  result  can  be  obtained. 

In  Section  II,  it  was  stated  that  the  function 
R(x)  describes  the  wearout  period  in  bath  tub  type  curves: 

R(*>  - rhn?  ' 8 > 0 

If  a new  parameter,  y , is  defined  which  is  between 
0 and  1 and  R’ (x)  is  now  defined  to  be  the  yth  power  of 
R(x)  : 
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R'  (x) 


(1  + 8x) 


6 > 0,  0 < y < 1 


(4.20) 


then  R' (x)  will  provide  for  a slowly  increasing,  rather 
than  rapidly  increasing,  in  wearout  period.  With  this 
revision  the  model  (2.6)  now  becomes 


z « G(x)  * xL(x)*R'(x) 


(1  + ax)  (1  + 8x) 


(4.21) 


(4.22) 


where  a > 0,  8 > 0 and  0 < y < 1. 

However,  after  we  put  another  variable  in  the  model, 
it  will  be  too  difficult  to  handle  in  the  previous  estimation 
technique  for  obtaining  values  of  the  parameters  a,  8,  A 
and  y because  of  nonlinearity  and  indefiniteness  of  the 
expression  (4.22)  . 

For  this  reason,  y will  be  assumed  constant  in  the 
previous  estimation  procedure.  Then  it  will  be  computa- 
tionally convenient.  Actually,  the  estimation  procedure 
does  not  require  any  change.  The  nonlinear  least  square 
estimation  approach  will  simply  be  used  for  different  values 
of  y.  Table  XIII  and  Fig.  15  demonstrate  the  parameters 
estimated  and  the  resulting  objective  function  values.  From 
Table  XIII,  it  can  be  easily  seen  that  the  new  result  has  a 
smaller  sum  of  squared  errors.  The  best  value  of  y seems 
to  be  near  0.95. 
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TABLE  XI 


R HUMAN  LIFE  BY  USING  SAME  INITIAL  VALUES 


TABLE  XIII 

ESTIMATION  VALUES  OF  a, 
USING  CONSTANT  y 


1.000 
.97 
0.95 
.93 
.90 
.85 
.7 


0.000834 

0.000512 

0.000363 

0.000740 

0.000685 

0.000592 

0.000319 


a 

6 

Objective 

Value 

0.42950 

0.01029 

4,335.89 

0.04051 

0.01182 

3,119.85 

0.01928 

0.01314 

2,494.10 

1.06371 

0.01352 

2,999.11 

1.81283 

0.01554 

3,267.47 

5.74169 

0.02031 

5,099.37 

3.17467 

0.04580 

1 

15,310.06 

C.  CONCLUSIONS 

Simple  analytical  hazard  models  have  been  developed 
and  fitted  to  situations  (data)  that  exhibit  bath  tub 
shaped  hazard  functions.  That  is,  failure  rates  may  be 
high  at  early  ages  ("infant  mortality") , constant  at 
intermediate  ages,  and  high  again  for  later  ages  ("wearout"). 
The  procedure  emphasizes  representations  of  the  inverse 
distribution  function;  simulation  is  thus  facilitated. 

The  failure  time  distributions  so  derived  should  be 
useful  in  analyzing  maintenance  and  replacement  policies. 

A least  squares  technique  for  fitting  the  hazard  models 
to  data  are  suggested  and  applied. 
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